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Abstract 

The critical properties of a cellular automaton model describing 
the spreading of infection of the Herpes Simplex Virus in corneal tissue 
are investigated through the dynamic Monte Carlo method. The model 
takes into account different cell susceptibilities to the viral infection, as 
suggested by experimental findings. In a two-dimensional square lattice 
the sites are associated to two distinct types of cells, namely, permissive 
and resistant to the infection. While a permissive cell becomes infected 
in the presence of a single infected cell in its neighborhood, a resistant 
cell needs to be surrounded by at least R > 1 infected or dead cells 
in order to become infected. The infection is followed by the death of 
the cells resulting in ulcers whose forms may be dendritic (self-limited 
clusters) or amoeboid (percolating clusters) depending on the degree of 
resistance R of the resistant cells as well as on the density of permissive 
cells in the healthy tissue. We show that a phase transition between 
these two regimes occurs only for i? > 5 and, in addition, that the 
phase-transition is in the universality class of the ordinary percolation. 



1 Introduction 



One of the most common and intensively studied diseases among humans is the 
herpes simplex virus (HSV) infection. Apparently, the unique symbiosis that exists 
in nature between humans and the HSV allows the viral particles to remain inactive 
(latent infection) in the cranial nerve ganglia after a primary infection, producing 
frequently recurring localized infections during the host's lifetime [0. The reacti- 
vation of HSV from latency may occur at any time and it is characterized by active 
viral replication in the epithelium causing vesicular eruptions in human mucosae 
and skin. The rupture of these vesicles and the consequent cell necrosis leave the 
characteristic herpetic lesion or ulcer. 

Basically, there are two distinct types of herpes simplex virus, namely, HSV 
Type I and HSV Type H. The former generally involves infection above the waist 
(ocular and facial) while the latter infects tissues below the waist. Here we discuss 
a mathematical model proposed to describe the growth of corneal ulcers caused by 
HSV Type I This infection is common and frequently causes corneal opacifi- 
cation. Traditionally the morphology of the corneal ulcers has been described as 
either dendritic or amoeboid. The dendritic ulcers are by far the more frequent 
form and, though they are self-limited in general, ocasionally they can enlarge pro- 
gressively changing to the amoeboid form. This is actually the natural course of the 
infection in the case of immunocompromised hosts or of inappropriate use of topic 
corticosteroids. In general, the amoeboid ulcers have a prolonged clinical course 
when compared to the dendritic ones. Regardless of their morphology, the ulcers 
are epithelial lesions that extend through the basement membrane whose swollen 
epithelial borders contain active viral particles. 

In order to carry out a more quantitative study of the ulcer morphology, the 
fractal dimension of clinically diagnosed HSV ulcers (including both dendritic and 
amoeboid forms) have been estimated suggesting that their outlines are fractal 
objects While the dendritic ulcers include branching and linear lesions, the 
geographic ulcers are no longer linear and as they increase in size, their perimeters 
become less and less irregular. In addition to its usefulness as a classification 
tool, the fractal properties of the ulcers may give information on the underlying 
mechanisms of viral spread within the epithelial tissue. For instance, a theory 
based only on the neurotropism of HSV and the dendritic-like distribution of nerve 
terminals can explain the branching pattern observed in dendritic ulcers Q , but it 
fails to explain the decrease of the fractal dimension (perimeter irregularity) with 
increasing ulcer sizes. An alternative explanation put forward by Landini 
which will be the main focus of this paper, considers the ulcer shape as the natural 
outcome of the contiguous spread of viral particles modulated by variations in the 
cell susceptibilities to infection. To take into account the fact that viruses only 
infect cells that have appropriate receptor molecules on their surface, those authors 
proposed a cellular automaton model for the HSV I spread in which the corneal 
epithelial tissue is modeled by a two-dimensional lattice. In their model, each 
lattice site may be occupied cither by a permissive cell (with probability q) or by 
a resistant cell (with probability 1 ~ q). More pointedly, a permissive cell becomes 
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infected whenever there are at least one infected cell in its neighborhood, while a 
resistant cell becomes infected if the number of infected and dead neighboring cells 
is larger than or equal to the integer parameter R > 1 that measures the degree of 
resistance of the cell ||^. 

The simulated ulcers obtained with the cellular automaton have the same qual- 
itative features of the clinical lesions and, in addition, for appropriate choices of 
the degree of resistance R a dramatic change on the morphology of the ulcers is 
observed as the initial concentration of permissive cells q increases beyond a cer- 
tain value This phenomenon was conjectured to be of a (qualitatively) similar 
nature as the ordinary percolation phase transition. The main contribution of this 
paper is to show, through the calculation of the dynamic and static critical expo- 
nents, that in the cases where a phase-transition does occur (i? > 5), the transition 
belongs indeed to the universality class of the ordinary percolation |6| . To carry 
out this analysis we use the so-called dynamic Monte Carlo method or spreading 
analysis [0, ^ whose idea is to study the spreading of the infection starting from 
a configuration with a single infected cell on the center of the lattice. Clearly this 
technique is very well suited to our investigation since the characterization of the 
spreading behavior of the infection is exactly the issue we address in this paper. 

The remainder of the paper is organized as follows. Following Landini ei a/ Q, 
m Sec. I we give the set of rules that govern the evolution of the HSV I infection 
in a two-dimensional square lattice and present the evidences for the existence of 
a threshold phenomenon or phase transition for R > 5. In Sec. || we characterize 
this phase transition using the dynamic Monte Carlo method which allows the 
computation of the critical dynamic exponents that describe quantitatively the 
spreading of the infection from a single infected cell. Finally, some concluding 
remarks are presented in Sec. 0. 



2 Model 

The cellular automaton model is defined in a square lattice consisting of {L + 
1) X (L -I- 1) sites, where each site is associated to a cell. Each cell is modeled by 
a four-state automaton corresponding to the different states of this cell: healthy 
permissive, healthy resistant, infected and dead. Except for the central cell, the 
initial state of any cell in the lattice is set either as permissive or resistant with 
probabilities q and 1 — q, respectively, so that there are no dead cells at the outset. 
The infection spreads from the single central infected cell and the ulcer (i.e., the 
cluster of dead cells) grows according to the following deterministic rules ||] : 

(1) An infected cell dies in the next time step. 

(2) A healthy permissive cell becomes infected if at least one of its neighboring 
cells is infected. 

(3) A healthy resistant cell becomes infected if at least i? > 1 of its neighboring 
cells are infected or dead. 
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Figure 1: Percolation probability 11 as a function of the initial density of 
permissive cells q for L = 1001 and (left to right) i? = 2, 3, 4, 5 and 6. 



The neighborhood of a given cell consists of its first and second nearest neighbors 
(Moore neighborhood). The infection and subsequent death of a resistant cell 
surrounded by R or more dead cells is justified by the lack of tissue support. In 
addition, this is necessary to prevent the occurrence of large ulcers with small islands 
of resistant cells, which are not observed clinically The four-state automaton 
considered allows transitions of one of the healthy states to the infected and dead 
states in a cyclic manner. At each time step we perform a parallel updating of all 
cell states. 

For R < 8 we are dealing with a variant of the so-called diffusion percolation 
process where the geometry changes via a dynamic process and the nature of the 
growth depends on the local environment For finite lattice sizes and open 
boundary conditions the above rules are repeated until either there are no more 
cells to infect or an infected cell reaches the lattice boundary. These different modes 
of termination generate dendritic (self-limited) and amoeboid (irrestricted) ulcers, 
respectively. It is interesting to note that the ordinary site percolation process is 
recovered for R > 8, since in this case a resistant cell can never become infected 
and so the infection can propagate only through the permissive cells. 
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To illustrate the dependence of the different termination modes, and hence ulcer 
forms, on the control parameters R and q of the model we present in Fig. ^ the 
fraction of irrestricted ulcers generated in 1000 runs. Each run corresponds to a 
different initial configuration of the lattice. Clearly, this fraction can be identified 
with the percolation probability 11 ||^, ^. Rather interestingly, we have found 
that the results for i? > 6 are indistinguishable within the numerical precision. 
Actually, this is expected since there is a preferred direction for the propagation 
of the infection, namely, from the center to the lattice boundaries, and so only the 
neighborhood facing the infection front matters to update the state of a healthy 
cell. Since the largest size of the error bars in this as well as in the next figure 
is twice the size of the symbols, they were ommited for the sake of clarity. More 
importantly, we have found that for i? < 4, the results become independent of the 
lattice size already for L > 101. 

However, for i? > 5 the dependence on the lattice size, illustrated in Fig. ^ for 
R — 5, indicates the occurrence of an atypical threshold phenomenon at a critical 
value Qc in the limit L oo. In fact, as q increases from to 1 the percolation 
probability 11 vanishes for q < q^, undergoes a discontinuous transition to some 
value H = He > at q ~ qc and then increases monotonically towards 1. This 
transition is atypical in the sense that He is not equal to 1 above qc, as in the case 
of the ordinary percolation transition ||^, ^ , which means that in this regime there 
is a finite probability that the infection does not percolate, i.e., a dendritic ulcer is 
formed. The reason for that is due simply to the fact that the spreading process 
starts from a single central cell so that if the infection happens to percolate in a 
lattice of a given size then it is certain to percolate in a smaller lattice too, i.e., 
n(Li) > U{L2) for Li < L2. In particular, n(3) = 1 - (1 -p)" yields an upper 
bound to n(oo). Of course, if the initial setting is such that there is an extensive 
number of infected cells, say aL with a < 1, randomly distributed over the bottom 
side of the lattice and periodic boundary conditions on the lateral sides, then the 
usual result 11^ = 1 is recovered |l^. In fact, since the curves for different lattice 
sizes do not cross, the standard finite size scaling analysis aiming at determining 
both and the spatial correlation length exponent v^^ for R > 5 (see, e.g., Ref. 
[^) fails spectacularly and so we have to resort to other means to estimate those 
quantities. 

3 Spreading analysis 

We turn now to the analysis of the spreading behavior of the viral infection starting 
from a single infected cell located in the center of a lattice of infinite size. Finite 
size effects are absent because the lattice size is taken large enough so that during 
the time we follow the evolution the infection front can never reach the lattice 
boundaries. This of course sets an upper limit to the time we can follow the viral 
spread and so, for instance, for lattices of size L = 4005 we let the infection evolve 
up to t = 2000. As usual, we concentrate on the time dependence of the following 
key quantities 0: (i) the average number of dead and infected cells n(t); (ii) the 
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Figure 2: Percolation probability 11 as a function of the initial density of 
permissive cells g for i? = 5 and L = lOl(O), 401(A), 701(v) and lOOl(x). 
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survival probability of the infection p{t) ; and (iii) the average mean-square distance 
over which the ulcer has spread r^it). For each time t we carry out 10'* independent 
runs, hence p{t) is simply the fraction of runs for which there is at least one infected 
cell in the lattice at time t. At the transition point qc we expect that the measured 
quantities obey the following scaling laws Q 

p{t) ~ t-' (1) 
n(t) - f (2) 
r^(t) ^ f- (3) 

where (5, t] and z are dynamic exponents. Since the fractal dimension d/ of the 
ulcer at a given time t is defined as nit) ^ r'^^ we have 

df = 2^ (4) 

at the critical point. Note that this equation is different from the one used in the 
studies of directed percolation (see, e.g., Ref. |ll)) because in the present case all 
runs generate an ulcer and so n{t) as well as r'^{t) are averages taken over all runs. 

In Figs. §,|4|and|| we present log- log plots oip{t), n{t) and r^(i), respectively, as 
functions of t in the vicinity of the critical point for R = 5. The asymptotic straight 
lines observed in these figures are the signature of critical behavior while upward 
and downward deviations indicate supercritical {q > qc) and subcritical {q < qc) 
behaviors, respectively. We recall that in the subcritical regime only dendritic ulcers 
are formed, while in the supercritical regime the formation of amoeboid ulcers is 
much more frequent (see Fig. |). The data shown in Fig. | yield = 0.3945±0.0002 
where the error is estimated by determining two values of q as close as possible to 
the critical point for which upward and downward deviations can be observed. A 
precise estimate for the dynamic critical exponents is obtained by considering the 
local slopes of the curves shown in the previous figures. For instance, the local 
slope S(t) is defined by |ll] 

\n[p{t)/p{t/8)] 
"^(') = h^8 ' 

which for large t behaves as 

d{t)^S+^ (6) 

where a is a constant. Analogous expressions hold for ri{t) and z{t). Hence plots 
of the local slopes as functions of 1/t allow the calculation of the critical ex- 
ponents. Applying this procedure for the critical curves we find the exponents 
S = 0.0870 ±0.0001, ri = 1.5866 ±0.0007 and z = 1.6843 ±0.0003. The errors in the 
critical exponents are, as usual, the statistical errors obtained by fitting the local 
slopes by straight lines in the large t regime. We expect, of course, that the (un- 
controlled) systematic errors are much larger than those. Using Eq. (^ we obtain 
df — 1.8840±0.0005 which is in very good agreement with the analytical prediction 
for the ordinary percolation df = 91/48 « 1.896 ^. 
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Figure 3: The log-log plot of p{t) as a function oi t ioi R = 5 and (top to 
bottom) q = 0.4, 0.395, 0.3947, 0.3945, 0.394, 0.393 and 0.39. 
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Figure 4: Same as fig. ^ but for n{t). 
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As the dynamic exponents S, rj and z for the ordinary percolation problem 
are not very well known, to show unambiguously that this ulcer formation model 
belongs to the universality class of the ordinary percolation we ought to estimate 
the static exponents f3 and i'±. We recall that the exponent (3 gives a measure 
of how the fraction of lattice cells belonging to an infinite cluster vanishes as the 
percolation threshold Qc is approached in the supercritical regime while is the 
correlation-length exponent in the space direction. To do so we calculate first the 
exponent which governs the decay of the concentration of infected cells i(t) in 
the subcritical regime. In fact, since in this regime the correlations are short-ranged 
one expects i{t) to decay exponentially 

l{t)^ Aiq)exp[-{q,-qpt] t ^ oo (7) 

where A (q) is some time independent function. Fig. ^ not only illustrates the 
adequacy of this assumption but permits also the evaluation of the decay constant 

A = (<7c - qP (8) 

from the asymptotic slopes of the curves Ini vs. t. The results presented in Fig. 0, 
showing the dependence of A on the distance qc — q from the critical point, allows 
the calculation of the exponent i^y as the slope of the straight line, yielding i^y = 
1.54 ± 0.03. Once this exponent is known we can use the scaling relations P — ly^^ 6 
and iy± — zv\\l2 to estimate the static exponents. We find /? ~ 0.134±0.003 and 

= 1.30 ±0.03 which, within error bars, are in agreement with the exact values of 
the corresponding exponents of the ordinary percolation, namely, /3 = 5/36 ~ 0.139 
and =4/3« 1.333 |, |. 

We have carried out a similar spreading analysis for i? > 6 and, as hinted in 
Fig. |l|, we have found a slightly larger percolation threshold, namely, qc = 0.4075 ± 
0.0002 which, within error bars, is shown to be independent of the value of i? > 6. 
Furthermore, since the larger the resistance parameter R, the more similar the 
ulcer formation problem is to the ordinary site percolation, we have found the same 
dynamic and static critical exponents as for the case i? = 5, as expected. 

4 Conclusion 

Using the dynamic Monte Carlo method we have shown unambiguously that the 
phase transition observed in the model for formation of herpes simplex ulcers pro- 
posed by Landini et al Q belongs to the universality class of the ordinary percola- 
tion. The value of this finding should not be underrated since the infection process 
actually resembles a diffusion percolation process where the growth depends on the 
local environment, in the sense that the decision on whether or not a resistant cell 
will become infected depends on the time-dependent states of several of its neigh- 
bors. Furthermore, since the ulcer formation model described here may be thought 
of as a damage spreading process, one could expect that the transition were in the 
universality class of the (2 -t- 1) directed percolation instead. However, as pointed 
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Figure 6: The log-linear plot of i(t) against t ioi R = 5 and (top to bottom) 
qc-q = 0.009, 0.011, 0.014, 0.016, 0.019, 0.021, 0.024, 0.027, 0.029 and 0.034. 
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Figure 7: The log-log plot of the time decay constant A against Qc — q for 
R = 5. The slope of the straight line yields !=a 1.54 zt 0.03. 
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out by Grassberger p2| , this is not so because in the ulcer formation model the 
damage never heals (even if it does not spread), i.e., the probability that an infected 
or dead cell becomes healthy is zero. 

The finding that for i? < 5 the model does not present a phase transition 
refiects the nontrivial role played by the resistance parameter R in this percolation 
process. In these noncritical cases the probability that an infinite or irrestricted 
ulcer is generated is given by the smooth size-independent curves shown in Fig. 
^ Similarly to a noncritical forest fire model the growing of this type of 

ulcer may be characterized by infection fronts with fractal dimension D whose value 
probably depends on the resistance parameter R (of course, D — \ ioY R— 1). An 
additional feature that makes the quantitative study of this viral spreading model 
rather challenging is the result that the percolation probability curves for different 
lattice sizes do not cross (see Fig. ||), which complicates enormously the estimate 
of the percolation threshold and critical exponents through the standard finite size 
scaling method. 

Some remarks on the biological interpretation of our results are in order. Ac- 
cording to the specialized literature ||, ^, ^, amoeboid ulcers are, in general, 
observed in immunocompromised patients or in patients that made inappropriated 
use of corticosteroids. In the present model, these conditions would correspond to 
a decrease of the degree of resistance R of the resistant cells or to an increase of 
the initial concentration q of permisssive cells. Although this model does not take 
into account the recurrent characteristic of this kind of infection, in which case 
the variability of q would probably play an important role, nor the possibility of 
variation of R during the course of the infection, its predictions are in qualitative 
agreement with the clinical observations. In fact. Fig. |^ points out the prevalence 
of amoeboid ulcers when R decreases or q increases. This agreement lends sup- 
port to the hypothesis that the morphology of the ulcers is determined by the viral 
spreading through cells with different susceptibilities to infection. 

To conclude we should mention that an extension of the original model pro- 
posed by Landini et al in which both the regeneration of dead cells as well as the 
spontaneous outbreak of infection anywhere in the lattice are taken into account 
has already been considered in the literature jl^. Interestingly, in this case the 
viral spreading model becomes very similar to the critical forest fire model with 
immune trees 0. In particular, the resistance parameter R of the ulcer for- 
mation model is akin to the immunity probability, i.e., the probability that a tree 
is not ignited though one of its neighbors is burning. According to a conjecture 
put forward by Grassberger [|l2j, the extended ulcer formation model should be in 
the universality class of directed percolation, since it allows for the regeneration of 
dead cells. This suggestion is strengthened by the finding that the forest fire model 
with immune trees is in that class of universality [ pTf . 
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